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ABSTRACT 

The X-ray spectra of accretion discs of eight stellar-mass black holes have been analyzed 
to date using the thermal continuum fitting method, and the spectral fits have been used 
to estimate the spin parameters of the black holes. However, the underlying model used 
in this method of estimating spin is the general relativistic thin-disc model of Novikov 
& Thorne, which is only valid for razor-thin discs. We therefore expect errors in the 
measured values of spin due to inadequacies in the theoretical model. We investigate 
this issue by computing spectra of numerically calculated models of thin accretion discs 
around black holes, obtained via three-dimensional general relativistic magnetohydrody- 
namic (GRMHD) simulations. We apply the continuum fitting method to these computed 
spectra to estimate the black hole spins and check how closely the values match the actual 
spin used in the GRMHD simulations. We find that the error in the dimensionless spin pa- 
rameter is up to about 0.2 for a non-spinning black hole, depending on the inclination. For 
black holes with spins of 0.7, 0.9 and 0.98, the errors are up to about 0.1, 0.03 and 0.01 
respectively. These errors are comparable to or smaller than those arising from current 
levels of observational uncertainty. Furthermore, we estimate that the GRMHD simulated 
discs from which these error estimates are obtained correspond to effective disc lumi- 
nosities of about 0.4 — 0.7 Eddington, and that the errors will be smaller for discs with 
luminosities of 0.3 Eddington or less, which are used in the continuum-fitting method. 
We thus conclude that use of the Novikov-Thorne thin-disc model does not presently limit 
the accuracy of the continuum-fitting method of measuring black hole spin. 

Key words: accretion, accretion discs; black hole physics; (magnetohydrodynamics) 
MHD; methods: numerical; relativity; X-rays: binaries 



1 INTRODUCTION 

Astrophysical black holes are described by just two parameters: 
their mass M and angular momentum J, with the latter usu- 
ally expressed in terms of the dimensionless spin parameter 
a t = cJ/GM 2 . While the mass M is relatively straightforward 
to obtain using dynamical measurements, the spin parameter a 4 
is less so. In accreting black holes, however, emission from the 
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inner disc gives us a handle on the spin. According to the model 
developed by Novikov & Thorne (1973, hereafter NT) 1 for a 
razor-thin accretion disc around a black hole, viscous evolution 
causes the accreting matter to move slowly inward along nearly 
Keplerian orbits until reaching the radius of the innermost sta- 
ble circular orbit (ISCO), after which the gas plunges into the 
black hole on a dynamical timescale. Thus, the inner edge of the 
viscous accretion disc is predicted to be very close to the ISCO. 
This link between the radius of the ISCO r ISC0 and the radius of 
the inner edge of the disc r in is well supported by empirical evi- 



1 This is the relativistic generalization of the standard thin-disc model 
of Shakura & Sunyaev (1973). 
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dence that the inner radius is constant in disc-dominated states 
of black hole binaries (McClintock et al. 2007; Steiner et al. 
2010a; and references therein), and by recent GRMHD simu- 
lations of thin accretion discs (Shafee et al. 2008a; Penna et al. 
2010; but see Noble et al. 2009, 2010). Therefore, measuring 
r ta gives one an estimate of r ISC0 . Since r ISC0 /M is a monotonic 
function of a„ (e.g., Shapiro & Teukolsky 1983), we can then 
obtain the spin of the black hole. This is the main technique 
currently being used to estimate the spins of stellar-mass black 
holes in binary systems. 

One of the ways of measuring r in involves 2 fitting the ther- 
mal X-ray continuum spectrum from the disc with the NT model 
spectrum (e.g., Zhang et al. 1997; Shafee et al. 2006; Davis et al. 
2006; Gou et al. 2009, 2010; Steiner et al. 2009, 2010a) using 
models such as KERRBB (Li et al. 2005) and BHSPEC (Davis 
& Hubeny 2006) in the data-analysis package XSPEC (Arnaud 
1996). From the fit, one obtains r to , or equivalently a,, if suit- 
able estimates of the black hole mass, inclination and distance 
are available (e.g., see Gou et al. 2009). Both KERRBB and BH- 
SPEC assume that the structure of the disc and its emission prop- 
erties are described accurately by the NT model. 

It is therefore clear that a crucial issue in black hole spin 
estimation via the continuum-fitting method is the NT model 
and its reliability. How much do real accretion discs with finite 
thickness differ from the NT disc? This question was addressed 
by Paczyriski (2000) and Afshordi & Paczyriski (2003), who ar- 
gued that deviations from the NT model decrease monotonically 
with decreasing disc thickness and that thin discs with dimen- 
sionless thickness \h/r \ <. 1 are well described by the model, if 
the viscosity parameter a<l. Their argument was confirmed 
by detailed calculations carried out by Shafee et al. (2008b). 
This still leaves open the question of whether magnetized discs 
might deviate substantially from the NT model even at small disc 
thicknesses (Krolik 1999). A number of recent studies of magne- 
tized discs using three-dimensional general relativistic magneto- 
hydrodynamic (3D GRMHD) simulations, including Shafee et al. 
(2008a), Noble et al. (2009, 2010) and Penna et al. (2010), have 
explored this question. These authors estimate that the luminos- 
ity and stress of the inner regions of simulated discs differ from 
the NT model by factors ranging from a few percent (Penna et al. 
2010) to as much as 20 percent (Noble et al. 2009). The ques- 
tion then is how much this departure affects measurements of 
black hole spin. 

We investigate this question using a very straightforward 
approach: we start with a disc model obtained via the above- 
mentioned GRMHD simulations (principally models similar to 
those described in Penna et al. 2010), compute the disc emission 
as a function of radius using a local blackbody approximation 
(assuming a constant spectral hardening factor), and use ray- 
tracing to compute the spectra. We then fit these spectra using 
KERRBB and compare the resulting spin estimate with the spin 
that was used in the GRMHD simulation. Our goal is similar to 
that of Shafee et al. (2008b), who performed the same analysis 
for a purely hydrodynamical disc, and of Li et al. (2010); the im- 
portant difference between the latter work and ours is that we 
use disc models obtained from GRMHD simulations, although 
we do not explore the effect of a finite photospheric height in 



2 Another method involves fitting the profile of the relativistically 
broadened iron line (e.g., Fabian et al. 1989, 2000; Laor 1991; Reynolds 
& Nowak 2003; Brenneman & Reynolds 2006). 



detail. Analogous work (though with a pseudo-Newtonian, not 
GRMHD, code) on the systematic errors in spin estimates ob- 
tained by fitting the broad iron emission lines from the inner 
accretion disc has been done by Reynolds & Fabian (2008) . 

We begin in §2 with a description of our method, and calcu- 
late in §3 the error in black hole spin estimates due to deviations 
of GRMHD discs from the NT model. We discuss in §4 the ob- 
servational uncertainties in black hole spin determination and 
compare these with the errors arising from use of the NT model. 
We conclude in §5 with a summary. Some technical details are 
discussed in Appendices A, B and C. 



2 METHOD 

2.1 Calculation of Disc Temperature and Velocity Profiles 
from GRMHD Simulations 

We work in Boyer-Lindquist (BL) coordinates x a = (t,r,6,cf)). 
We use geometric units where the speed of light c, the gravita- 
tional constant G and the Planck constant are set to unity, and 
measure all lengths and times in units of the black hole mass M. 

For this project we reran the three-dimensional general 
relativistic magnetohydrodynamic (3D GRMHD) simulations of 
thin discs described in Penna et al. (2010), for four values of 
the black hole spin: a t = 0, 0.7, 0.9, 0.98. For completeness, we 
briefly review the simulation setup here. The simulations solve 
the GRMHD equations for plasma around a rotating black hole 
using the code HARM (Gammie et al. 2003) with numerous re- 
cent improvements, including 3D capabilities (McKinney 2006; 
McKinney & Blandford 2009). The gas is initially in a torus in hy- 
drodynamic equilibrium surrounding the black hole (De Villiers 
et al. 2003; Gammie et al. 2003). The spin axes of the torus and 
the black hole are aligned. The torus is seeded with a magnetic 
field consisting of four poloidal loops arranged in the radial di- 
rection. We use a polytropic equation of state for the gas, p oc p r , 
where p, p and y are the pressure, density and adiabatic index 
respectively, and choose y = 4/3, as appropriate for a radiation 
pressure dominated disc. 

To keep the disc thin, we use a simple cooling prescrip- 
tion that drives the gas to its initial entropy on a dynamical 
timescale 3 . The energy removed by the cooling prescription is 
assumed to be completely lost by the accretion flow; it has no 
dynamical effect on the accreting gas (the energy lost to cooling 
is of course tracked and is later used to compute the disc lumi- 
nosity profiles shown in Fig. 1). The disc thickness is specified 
by the quantity \h/r\, where h is the density scale height of the 
disc above the midplane, \h\ = j p \z\ dz/ j p dz, and r is the 
cylindrical radius. Our simulated thin discs have \h/r\ = 0.05, 
0.04, 0.05 and 0.08 respectively for a„ = 0, 0.7, 0.9 and 0.98. 
Following Penna et al. (2010), we perform a temporal and az- 
imuthal average over the steady-state portion of the simulation 
results to average over the fluctuations introduced by turbu- 
lence, since we are interested in the mean behaviour of the ac- 
cretion flow. Finally, since our discs are geometrically thin, we 



3 One change from Penna et al. (2010) is that in the present work we 
cool all the gas, including the gas in the corona, whereas in most of their 
simulations, Penna et al. cooled only the disc region of the flow. The 
present simulations are similar to the "no-tapering model" described in 
§5.7 and Fig. 13 of their paper. 
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perform a density-weighted average in the polar direction to ob- 
tain the vertically-integrated disc structure. This process of col- 
lapsing the simulated disc into the equatorial plane circumvents 
the difficulty of defining a photosphere for the disc and calculat- 
ing emission profiles along it. A proper treatment would require 
a full radiative transfer calculation (e.g., Davis et al. 2005; S^- 
dowski et al. 2010), which is beyond the scope of this work. 

At the end of this process, both components required to 
calculate the spectra are available: the radial profile of the fluid 
four-velocity u M (r) in the equatorial plane and, from the energy 
removed by the cooling prescription 4 , the profile of the emit- 
ted flux F(r) = dE/{rdrd(f)dt) (where E is the energy emit- 
ted from one side of the disc as measured by an observer at 
infinity). When calculating the spectrum, the effect of electron 
scattering in the disc is taken into account indirectly using a 
colour correction (or spectral hardening) factor / col . For the re- 
sults presented here, we choose for simplicity a fiducial value of 
/ col = 1.7 (Shimura & Takahara 1995). Detailed models of disc 
atmospheres by Davis et al. (2005) and Davis & Hubeny (2006) 
indicate that / col can vary between 1.4 and 1.7, but this extra 
sophistication is not necessary for the simple tests described in 
the present paper. 

The flux profiles obtained from the GRMHD simulations are 
only reliable within a radius inside which the accretion flow has 
reached steady state. Outside this radius, which we call the in- 
flow equilibrium radius r ie (see Penna et al. 2010 for a defini- 
tion), we extend the profiles using the analytical disc model of 
Page & Thorne (1974). The procedure we use is described in 
Appendix A. Within a certain range, the exact choice of r ie does 
not affect the results of the extrapolation, as we show in that 
appendix. 

Fig. 1 compares the luminosity profiles d(L/M)/d(lnr) 
that we obtain to those in the standard NT disc model, for four 
values of the spin: a s = 0, 0.7, 0.9, 0.98. Here, M is the ac- 
cretion rate, and the luminosity I = 2dE/dt = 2 j Frdrdcj) = 
An j Frdr, so d(_L/M)/d (In r) = 4nr 2 F(r)/M (the extra factor 
of 2 is to account for emission from both sides of the disc) . 

The NT model has no radiation from inside the ISCO, 
whereas the simulations show some emission from this region. 
In addition, the peak of the emission in the simulated discs is 
seen to shift inward relative to the NT model. These effects are 
similar to those obtained by Sadowski (2009) for slim discs. As 
explained in that work, for large enough accretion rates (> 0.3 
Eddington), the accretion flow starts becoming radiatively inef- 
ficient at moderate radii (r ~ 10M — 30M); as a result, some 
of the heat generated by viscous dissipation at larger radii is ad- 
verted inward and released at smaller radii. Another important 
effect is that discs with finite thickness have a non-vanishing 
stress at the ISCO (in contrast to the razor-thin discs which the 
NT model considers for which the stress is expected to van- 
ish). This stress leads to additional viscous dissipation at radii 
r ~ r Isco . In our model, the inward shift in the emission peak 
due to both of these effects mimics a decrease in r ISC0 (see 
Fig. 1), i.e., an increase in the predicted black hole spin. As a 
result, fitting the GRMHD disc spectrum using the NT model 
leads to an overestimate of the black hole spin, as we shall see 
in §3. 



4 We only include the energy removed from the bound gas, since in- 
cluding all the gas results in an overestimate of the luminosity (Penna 
etal. 2010). 
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Figure 1. Luminosity profiles from the GRMHD simulations (solid lines) 
compared with those from the NT model (dashed lines) for a„ = 0, 0.7, 
0.9 and 0.98 (bottom to top). The disc thicknesses are \h/r\ = 0.05, 
0.04, 0.05 and 0.08 respectively for these runs. The ISCO is located at 
the radius where the NT disc luminosity goes to zero. 



2.2 Calculation of the Spectra 

To calculate the spectrum, we assume that the flux F(r) is emit- 
ted in the form of colour-corrected blackbody radiation (/ col = 
1.7), either isotropically or with limb-darkening, as seen in the 
comoving frame of the fluid. We use a standard limb-darkening 
prescription (Eq. 5 below). We assume that after emission, the 
radiation propagates in vacuum. 

Were the accretion disc non-relativistic, the calculation of 
the spectrum would be almost trivial (see, e.g., Frank et al. 
2002); one would divide the disc into annuli, define an effec- 
tive blackbody temperature T eff (r) = [F(r)/cr] 1/4 in each annu- 
lus, where a is the Stefan-Boltzmann constant, use the tempera- 
ture and colour-correction factor to obtain the specific intensity 
J v disc (r) of the emitted radiation at the disc surface, and inte- 
grate it over the disc surface to obtain the observed spectrum. 

Relativity introduces three complications: (1) The effective 
temperature has to be defined in the comoving frame of the 
fluid, and so we need to transform F(r) from the BL frame into 
the comoving frame. (2) Redshift between the comoving frame 
and the observer's frame, both gravitational and due to Doppler 
boosting, has to be taken into account. Since the photon paths 
around a black hole are complicated, the direction in which a ray 
needs to be emitted in the comoving frame such that it reaches 
the observer is not known a priori, which is a problem for the 
redshift calculation. (3) One needs to know the emission direc- 
tion to take limb darkening into account as well. Points (2) and 
(3) require integrating the geodesic equations to calculate the 
photon paths, which is usually referred to as "ray-tracing." This 
approach has been applied extensively in the literature to a va- 
riety of problems, starting with Cunningham & Bardeen (1973) 
and Cunningham (1975) (see Dexter & Agol 2009 and refer- 
ences therein). In particular, KERRBB (Li et al. 2005) uses this 
technique to compute thin-disc spectra. 

We perform ray-tracing numerically using the routines 
developed by Shcherbakov & Huang (2011) and applied in 
Shcherbakov et al. (2010). We choose a line of sight to the ob- 
server with an inclination angle of i relative to the black hole 
spin axis. At a sufficiently large distance from the black hole 
(r ~ 10 5 ) we set up an image plane perpendicular to the line 
of sight and shoot rays from it parallel to the line of sight. We 
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Figure 2. Spectra from the simulated (solid line) and NT (dashed line) 
discs, for a t = 0.9 and i = 75°. 



follow these rays until they hit the disc 5 , by directly integrating 
the (second-order) geodesic equations: 



d 2 x a ^ dx p dx r 
dX 2 +T ^~dX~dX 



: 0, 



(1) 



where A is an affine parameter along the geodesic, and r| are 
the connection coefficients. The aim is to obtain the specific in- 
tensity of each ray, which can then be integrated over the 
image plane to obtain the observed flux: 



D 2 



LdA. 



(2) 



Since I v /v 3 is a Lorentz invariant, we can immediately relate 
the specific intensity 7 r in the image plane to the intensity I v com 
in the comoving frame of the fluid at the point of emission: 

\3 . 



A', com (^/^com) — ^v,comX 

where % = v/v com is the redshift factor, and 



2/", 4 v 3 

■> col COI 



-T. 



(3) 



(4) 



v * exp(v com /fc B / co iT com )- 

Here, / col is the spectral hardening factor mentioned earlier, 
which we set to 1.7 in this work, and T is the limb-darkening 
factor (see, e.g., Li et al. 2005): 



T = { 



So finally we have 
L = 



2/ c > 3 



exp(v/fc B / col ^T com )- 1 



isotropic emission 
limb-darkened emission 



T, 



(5) 
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Figure 3. High-energy portion of the spectra from the simulated (solid 
lines) and NT (dashed lines) discs, for i = 75° and three values of the 
black hole spin: a, = 0, 0.7, 0.9. 



Thus, to calculate the spectrum, we need the effective temper- 
ature in the comoving frame r com , the redshift factor % an d the 
angle 6 com between the emitted ray and the disc normal in the 
comoving frame. The first of these is obtained by transforming 
the emitted flux F(r), which is initially calculated in the Boyer- 
Lindquist frame, into the comoving frame, and the last two by 
transforming the ray four-momentum. We show the details in 
Appendix B. 

To calculate spectra using Eq. 7, we use the following fidu- 
cial parameters: black hole mass M = 10M S , accretion rate 
M = 0.1M Edd , and distance to the black hole D = 10 kpc. We 
choose a spectral energy range of 0.1 keV to 10 key divided into 
1000 logarithmically spaced bins. Fig. 2 compares the spectra 
from the simulated and NT discs for a t = 0.9 and i = 75°. The 
peak of the spectrum of the simulated disc is shifted to a slightly 
higher energy relative to the NT spectrum, and the peak flux is 
also higher. This is precisely the effect that increasing the black 
hole spin would have on the NT spectrum, and as we shall see 
in §3, fitting the simulated spectra leads one to overestimate the 
spin. 

Fig. 3 and Fig. 4 show what happens to the difference be- 
tween the simulated and NT spectra when we vary the spin 
and the inclination respectively. The effect is visible at the high- 
energy end of the spectrum. The effect of the inclination is very 
strong. This is because of the excess luminosity from the inner 
region of the simulated disc; this excess is more noticeable at 
higher inclinations due to beaming of the emitted radiation. 



2/ c > 3 



exp(v/fc B / col ^T com )-l 



-TdA 



(7) 



5 This is more straightforward than shooting rays from the disc, since 
as mentioned earlier, the direction in which the rays need to be emitted 
from the disc such that they reach the observer is not known a priori. 
This approach was pioneered by Marck (1996); see also Hameury et al. 
(1994). 



2.3 Tests 

We tested our code by comparing the spectra it produces for 
an NT disc to those produced by KERRBB itself using the 
f akeit command in XSPEC, for the following range of pa- 
rameters: black hole masses of 5, 10, 15 M , spin parameters 
a t = 0,0.7,0.9, observer inclinations i = 15°, 45°, 75°, accre- 
tion rates M/M Edd = 0.1,0.2 and distances D =10 kpc, 20 kpc, 
with and without using limb darkening. At a grid resolution of 
N b x N p = 100 x 100 (see the description of our grid in Ap- 



© 0000 RAS, MNRAS 000, 000-000 



BH Spin from Realistic Temperature Profiles 5 



! = 15° 

i = 15°, NT 



i = 45° 
i = 45°, NT 



i = 75° 
i = 75°, NT 



3 

E (keV) 



10 



Figure 4. High-energy portion of the spectra from the simulated (solid 
lines) and NT (dashed lines) discs, for a t = 0.9 and three inclinations: 
i = 15°, 45°, 75°. 



Table 1. Spin estimates obtained by fitting the simulated spectra (for 
limb-darkened emission) with RERRBB, for a range of spins a, and ob- 
server inclination angles i. The model identified as "1 loop" corresponds 
to a GRMHD simulation that has one poloidal magnetic loop in its ini- 
tial disc configuration; all the other models have four loops arranged 
radially. 





a, = 


a, = 0, 1 loop 


a, =0.7 




\h/r\ =0.05 


\h/r\ = 0.07 


|h/r| = 0.04 


i = 0° 


0.08 ± 0.02 


0.06 ± 0.01 


0.71 ± 0.01 


£ = 15° 


0.08 ± 0.02 


0.06 ± 0.01 


0.72 ± 0.01 


i = 30° 


0.09 ± 0.02 


0.07 ± 0.02 


0.72 ± 0.01 


i = 45° 


0.10 ±0.02 


0.09 ± 0.02 


0.73 ± 0.01 


i = 60° 


0.11 ± 0.02 


0.18 ± 0.01 


0.76 ± 0.01 


i = 75° 


0.15 ± 0.04 


0.37 ± 0.01 


0.80 ± 0.02 




a, 


= 0.9 a t -- 


= 0.98 




\h/r\ 


= 0.05 \h/r 


= 0.08 



= 0° 

: 15° 

= 30° 
= 45° 
= 60° 
= 75° 



0.905 ± 0.002 
0.906 ± 0.002 
0.907 ± 0.003 
0.908 ± 0.003 
0.914 ± 0.005 
0.929 ± 0.006 



0.985 ± 0.001 
0.985 ± 0.001 
0.985 ± 0.001 
0.986 ± 0.001 
0.987 ± 0.001 
0.991 ± 0.001 



pendix C), the spectra calculated with our code converge to the 
KERRBB spectra in all these cases. This confirms that the code is 
robust. 



Table 2. Absolute and fractional errors in the estimated radius of the in- 
nermost stable circular orbit, r ISC0 , corresponding to the spin estimates 
in Table 1 . 



3 RESULTS 

We use the following fiducial parameters for our spectra, as 
mentioned earlier: black hole mass M = 1OM , accretion rate 
M = 0.1M Edd , and distance to the black hole D = 10 kpc. The 
spectra are fitted using KERRBB (without using returning radia- 
tion, since we do not include it in our spectrum calculation) to 
obtain the black hole spin a„ and accretion rate M. 

One more thing needs to be taken care of. Even though 
we average the GRMHD simulation results azimuthally and over 
time to remove the effects of turbulence and obtain a mean pro- 
file for the flux F(r) and the gas four-velocity u M (r), there is 
still some stochastic variation in the spin estimates with time. 
Therefore, we divide the steady-state portion of each simulation 
into chunks of duration At = 1000M (each simulation has 4-5 
such chunks), obtain a spin estimate from each chunk, and then 
quote the mean spin estimate and the error in the mean for each 
simulation 6 . 

The results for limb-darkened emission 7 are shown in 
columns 1, 3, 4 and 5 of Table 1. As expected, the fitted values 



In addition, there are a couple of potential sources of systematic er- 
ror: (i) our somewhat arbitrary choice of the matching radius used for 
extending the luminosity profiles beyond the inflow equilibrium radius 
(see Appendix A), and (ii) the fact that we restrict ourselves to the bound 
gas when calculating the luminosity profiles, as mentioned in footnote 4. 
We estimate the systematic error due to these two factors and, to be 
conservative, include them in quadrature in the error estimates that we 
quote in Tables 1 and 3. 

7 We also looked at spectra generated using isotropic emission. To fit 
these spectra we turned off the limb-darkening flag of KERRBB. The 
resulting spin estimates are very similar to those obtained using limb- 
darkened emission, so we do not show them here. 





a, =0 


a, = 0.7 




r ISCO = ^ 


r isco = 3-39 


! =0° 


-0.26 (-4.3%) 


-0.07 (-2.0%) 


! = 15° 


-0.27 (-4.5%) 


-0.07 (-2.2%) 


i = 30° 


-0.29 (-4.9%) 


-0.10 (-2.9%) 


i = 45° 


-0.33 (-5.4%) 


-0.15 (-4.5%) 


i = 60° 


-0.37 (-6.2%) 


-0.28 (-8.2%) 


i = 75° 


-0.49 (-8.2%) 


-0.51 (-15.1%) 




a, = 0.9 


a„ = 0.98 




r isco = 2-32 


risco = 1-61 


! =0° 


-0.04 (-1.6%) 


-0.07 (-4.1%) 


! = 15° 


-0.04 (-1.7%) 


-0.07 (-4.2%) 


i = 30° 


-0.04 (-1.9%) 


-0.07 (-4.5%) 


i = 45° 


-0.06 (-2.5%) 


-0.08 (-5.2%) 


i = 60° 


-0.10 (-4.1%) 


-0.11 (-6.9%) 


i = 75° 


-0.21 (-9.0%) 


-0.18 (-11.0%) 



are different than the ones used in the GRMHD simulations. The 
differences are largest at low spins. It is easy to understand why 
the difference is not constant; the dependence of the disc tem- 
perature profile (which determines the shape of the spectrum) 
on the spin is highly nonlinear. In particular, the position of the 
spectral peak in the NT model strongly depends on the radius 
of the ISCO (r ISC0 ), to the extent that one can think of KERRBB 
as fitting for r ISC0 instead of a„. There is a one-to-one relation- 
ship between a, and r ISC0 (see., e.g., Shapiro & Teukolsky 1983), 
shown in Fig. 5. At high spins, r ISC0 varies very rapidly as a func- 
tion of a„; conversely, a„ varies relatively slowly as a function 
of 

r isco- Thus, a given fractional error in r ISCO translates into a 
much smaller error in a„ at higher spins than at smaller ones. 
This is illustrated by the two gray bands in Fig. 5. Each band 
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Table 3. Fitted accretion rates for the cases in Table I, in units of 
10 18 gs -1 . The input values in the GRMHD simulations are denoted 
by M input . 



Figure 5. Relation between the black hole spin parameter a„ and the 
radius of the innermost stable circular orbit r Isco . 



represents a range of ±10% in r ISC0 , around r ISCO = 6 (upper 
band, corresponding to a t = 0) and r ISC0 = 1.61 (lower band, 
corresponding to a, = 0.98). The corresponding range in a„ is 
±0.2 at a» = 0, but only ±0.01 at a, = 0.98. 

For reference, we show the errors in the estimated r ISC0 
in Table 2 (the exact values of r ISCO for a t = 0, 0.7, 0.9 and 
0.98 are 6, 3.39, 2.32 and 1.61), and the fitted accretion rates 
in Table 3 (the fiducial rate of 0.1M Edd used here corresponds to 
(2.45, 1.35, 0.898, 0.5982) x 10 18 gs" 1 for the four spins, assum- 
ing that the accretion efficiencies are given by their NT values) . 
It is interesting to note that although the errors in the fitted spins 
are relatively large, the errors in the accretion rates are only a 
few percent. 

The other important effect is that of the observer inclina- 
tion: at high inclination, the error in the spin estimate is larger. 
This is because the difference between the disc temperature pro- 
files in the NT model and the simulations is significant only in 
the inner disc. At low inclination angles, the combined effect of 
gravitational redshift and beaming of the radiation (the latter of 
which concentrates the radiation close to the equatorial plane) 
results in this difference not being noticeable in the spectrum. At 
high inclination angles, on the contrary, beaming enhances the 
difference, causing the error in the fitted spin to increase. 

Changing the black hole mass, accretion rate or distance 
only changes the overall scaling of the spectrum; therefore, 
there is no effect on the shape of the spectrum or the spin esti- 
mates. We should note, however, that since the GRMHD simula- 
tions use dimensionless quantities, the gas mass scale in the sim- 
ulations is arbitrary. Therefore, we have the ability to choose any 
accretion rate for a given disc thickness. For real discs, this is cer- 
tainly not the case. The relation between the disc thickness and 
the luminosity is discussed in more detail in §3.1. We show there 
that the disc thicknesses used in our simulations Qh/r\ = 0.05, 
0.04, 0.05 and 0.08 for a t = 0, 0.7, 0.9 and 0.98 respectively) 
correspond to L/L Edd = 0.5, 0.4, 0.5 and 0.7 respectively; there- 
fore, strictly speaking, our estimates of the errors in the spin 
determination are only applicable for these luminosities. 

We carried out a test run at a t = with a different initial 
magnetic field configuration that has one poloidal loop instead 
of four as in our other runs. This model is closer in spirit to the 
simulations run by Noble et al. (2009, 2010). We find that the 
one-loop model gives hotter spectra and a larger error in the 
derived value of the spin at large inclination angles (compare 





a, = 


a* = 0.7 




M ta p Ut = 2.45 


Minput = 1-35 


i = 0° 


? 46 ± 01 


1 37 ± 01 


i = 15° 


J 46 ± 01 


1 37 ± 01 


i = 30° 


2 46 ± 01 


1 37 ± 01 


i = 45° 


2.47 ± 0.01 


1.36 ± 0.01 


i = 60° 


2.48 ± 0.01 


1.34 ± 0.02 


i = 75° 


2.49 ± 0.04 


1.30 ± 0.04 




a* = 0.9 


a* = 0.98 




Minput = 0.898 


M^ut = 0.5982 


i = 0° 


0.901 ± 0.001 


0.5983 ± 0.0003 


i = 15° 


0.901 ± 0.001 


0.5982 ± 0.0004 


i = 30° 


0.902 ± 0.001 


0.5983 ± 0.0005 


i = 45° 


0.903 ± 0.001 


0.5985 ± 0.0005 


i = 60° 


0.904 ± 0.003 


0.5980 ± 0.0008 


i = 75° 


0.892 ± 0.007 


0.594 ± 0.001 



the first two columns of Table 1). This agrees with the results 
described by Noble et al. (2010) and Penna et al. (2010), who 
investigated the behavior of other diagnostics such as the angu- 
lar momentum and shear stress and showed that GRMHD discs 
calculated from single-loop initial conditions generally deviate 
more strongly from the NT model compared to discs obtained 
from multi-loop initial conditions. Penna et al. (2010) argued 
that the multi-loop case is more natural since it better mimics 
disc turbulence, whereas the one-loop case might introduce an 
artificial long-range radial coherence in the solution. 

The errors in the spin estimates could be due to a number 
of reasons: (1) The disc emissivity profile outside the ISCO is 
different in the simulations compared to the NT model, as Fig. 1 
shows; (2) The simulations have some radiation coming from 
the plunging region inside the ISCO; and (3) Even outside the 
plunging region, the radial component of the gas four-velocity 
is not negligible. To find out which of these is the dominant 
effect, we calculated some spectra from the simulated discs by 
excluding the region inside the ISCO and setting the gas velocity 
outside the ISCO to its NT value. Any residual differences in the 
spin estimates would solely be due to (1). 

The results are shown in Table 4. We see that spin estimates 
obtained from the GRMHD simulations are still significantly dif- 
ferent from the true values. This shows that the dominant reason 
for the errors in the spin estimates is the fact that, even outside 
the ISCO, the disc emissivity profile in the simulations is differ- 
ent from the NT profile; more specifically, that the peak of the 
profile is shifted to smaller radii, as mentioned in §2. We should 
note, however, that for discs thicker than those considered in 
this work by about a factor of 2 or more, the effect of the plung- 
ing region is important. This follows from the finding of Penna 
et al. (2010) that deviations of the GRMHD simulations from 
the NT model increase with increasing disc thickness. 



3. 1 Effective Accretion Rates of the GRMHD Models 

The GRMHD models that we have used in this study make use of 
dimensionless quantities and do not include detailed radiation 
transfer. Hence there is no direct way of estimating the physi- 
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Table 4. Spin estimates from spectra obtained by excluding the plunging 
region and setting the four-velocity in the disc to its NT value (second 
column) compared with the original spin estimates from Table 1 (third 
column) . 



a t = 0.9, i = 75° 
a„ = 0.9, i = 45° 
a = 0, i = 75° 



0.92 0.93 
0.91 0.91 
0.13 0.15 



cal mass accretion rate (gs _1 ) or the true radiative luminosity 
(ergs -1 ) of the models. To estimate these quantities we use an 
indirect method, in which we compare the vertical thicknesses 
of the simulated discs against physical disc models that do in- 
clude radiation transfer and radiation pressure and solve for the 
vertical disc structure. 

We use two models for this comparison. One is a semi- 
analytical model of a slim disc (Sadowski et al. 2010) which 
goes beyond the NT model by including the effect of energy 
advection in the radial equations. At each radius r, the model 
solves the condition of vertical hydrostatic equilibrium and in- 
cludes radiative transfer approximately. The other model (Davis 
et al. 2005) assumes the NT model for the radial structure but 
carries out a careful and detailed computation of radiation trans- 
fer, including non-LTE effects, at each r. This model is identical 
to the XSPEC model BHSPEC (Davis & Hubeny 2006). Each of 
these models treats some part of the physics very well, but nei- 
ther has all the ingredients one would like to include, viz., ad- 
vection, full radiative transfer, magnetic fields, deviations from 
hydrostatic equilibrium, self-irradiation, etc. 

Fig. 6 shows the disc thickness as a function of r for a, = 0, 
as predicted by the slim-disc model (solid lines, correspond- 
ing to I/L Edd = 0.3, 0.4, 0.5, bottom to top) and BHSPEC 
(dashed lines, same set of luminosities). These curves should 
be compared with the disc thickness in the GRMHD simula- 
tion (dotted lines) . The top panel shows the density scale height 
\h\ = j p\z\dz/ j p dz, while the bottom panel shows the rms 

height h^^i jpz^z/jpdz) 1 ' 2 . 

At radii close to the ISCO, the slim disc and BHSPEC models 
indicate that the disc thickness plunges to small values whereas 
the GRMHD simulation shows a much smaller decrease. We be- 
lieve there are at least three reasons for this discrepancy: (i) the 
GRMHD simulations cool the gas by forcing it towards a con- 
stant entropy, which may not be justified in the plunging re- 
gion; (ii) the simulated GRMHD disc includes magnetic fields 
whose pressure provides additional support in the vertical di- 
rection whereas the other two models do not; and (iii) the sim- 
ulated disc begins to deviate from hydrostatic equilibrium as the 
radial velocity becomes large near the ISCO and the gas has less 
time to reach equilibrium, whereas the other models hardwire 
the condition of hydrostatic equilibrium at all radii. We estimate 
that the last two are only important well inside the ISCO, how- 
ever, while the discrepancy sets in already at larger radii. These 
are interesting issues which we hope to explore in the future. 
For the purposes of this section, we simply ignore the region of 
the simulation near the ISCO. 

For the comparisons described here, we select a radius of 
r = 12M = 2r ISC0 , which is well outside the ISCO, and deter- 
mine the luminosities at which the slim disc and BHSPEC models 
give the same disc thickness as we obtain in the simulated disc. 
We see from Fig. 6 that the thickness measure \h/r\ ~ 0.05 in 
the simulated GRMHD disc corresponds to L/L Eii ~ 0.5 accord- 




Figure 6. Profiles of \h/r\ for a t = from the slim-disc model of Sa- 
dowski et al. (2010) (solid lines) and the detailed radiative transfer 
model BHSPEC of Davis & Hubeny (2006) (dashed lines) for luminosi- 
ties L/iEdd = 0-3, 0.4, 0.5 (bottom to top in each panel), compared with 
the profile obtained from the GRMHD simulation (dotted lines) . The top 
panel shows the density scale height \h\ = J p |z| dz/ J p dz, while the 
bottom panel shows the rms height h lms = {J p z 2 dz / J p dz) 1 ' 2 . 



ing to the slim-disc model and ~ 0.4 according to BHSPEC. A 
comparison of the thickness measure h lms /r gives slightly larger 
values of L/L Mi . Similar analysis (again at r = 2r ISC0 ) for a„ = 
0.7, 0.9 and 0.98 shows that L/L Edd ~ 0.4, 0.5 and 0.7 respec- 
tively according to the slim-disc model, and 0.4, 0.4 and 0.6 
respectively according to BHSPEC. 

The above luminosities are much higher than the typical 
luminosities (L/L Edd < 0.3) at which black hole binary spectra 
are analyzed for spin determination using the continuum-fitting 
method (from Fig. 6, this corresponds to \h/r\ < 0.03). There- 
fore, the errors in the spin determination quoted in §3 are not 
directly applicable to observations. At lower luminosities, the 
accretion disc would be thinner, as Fig. 6 shows. Due to compu- 
tational resource requirements, it is not currently possible for us 
to perform GRMHD simulations for discs thinner than those pre- 
sented here. What we can do instead is to scale the results from 
the simulations to more realistic thinner discs at lower luminosi- 
ties. Table 5 compares the spins obtained from two GRMHD sim- 
ulations corresponding to a„ = 0, one with \h/r\ = 0.05 (which 
we have focused on so far) and another with \h/r\ = 0.09. We 
see that the error in the spin estimate is much larger for the lat- 
ter model. Thus, it is clear that at the lower luminosities that 
are interesting from an observational point of view, the errors in 
the spin estimates would be significantly smaller than those in 
Table 1. 

As an aside, we note that McClintock et al. (2006) calcu- 
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Table 5. Spin estimates for a„ = from GRMHD simulations of a thicker 
disc Qh/r\ = 0.09) compared with those for the thinner Qh/r\ = 0.05) 
disc shown in Table 1. 



\h/r\ 


i/^Edd 


i = 15° 


i = 45° 


i = 75° 


0.05 


0.5 


0.08 


0.10 


0.15 


0.09 


0.9 


0.08 


0.19 


0.50 



Table 6. Spin estimates from spectra calculated using a finite photo- 
spheric height for the disc (|h p h t/ r ! = 0-2) for both NT and simulated 
disc temperature profiles, at selected values of spin and inclination. 







NT disc 


simulated 
disc 


simulated disc 
(original estimates 
from Table 1) 


a, = 0.9, 


i = 75° 


0.88 


0.90 


0.93 


a, = 0.9, 


i = 45° 


0.89 


0.90 


0.91 


a, = 0, i 


= 75° 


-0.05 


0.09 


0.15 



lated the height of the disc photosphere as a function of radius 
for an NT disc (see Figure 17 in their paper). The disc heights 
they quote are larger than the density scale heights shown in 
Fig. 6 by roughly a factor of 2-3. We have calculated the lo- 
cation of the photosphere corresponding to the slim disc and 
BHSPEC models, and find that they agree fairly well with the 
values obtained by McClintock et al. (2006). 



3.2 Effect of a Finite Photospheric Height 

So far we have been calculating spectra using equatorial pro- 
files of the emitted flux F(r) and the fluid four velocity u M (r), 
which are obtained, as mentioned in §2, by vertically integrating 
the GRMHD simulated disc structure. This integration effectively 
collapses the disc into the equatorial plane, which is where the 
disc emission is assumed to originate from. The errors in the 
spin estimates quoted above have therefore been purely due to 
the departure of the equatorial flux and fluid velocity from their 
NT values. In reality, however, the observed disc emission comes 
from the photosphere, which is at a finite height above the equa- 
torial plane. The effect of this on the spin estimates needs to be 
checked. 

The method we use for doing this is very crude; our only 
goal here is to find out whether or not this effect could be impor- 
tant. The slim disc and BHSPEC models mentioned above, and 
the calculations of McClintock et al. (2006), show that the pho- 
tosphere height is about 2-3 scale heights. Since we want to see 
how large the effect of off-midplane emission could possibly be, 
we choose a larger photosphere height: h phot /r ~ 4|fr/r|, where 
\h/r\ is the disc half- thickness measured at one scale height 
above the disc midplane. We then repeat our ray-tracing compu- 
tation, but following the geodesies until they hit the photosphere 
instead of the equatorial plane (the photosphere as defined 
above corresponds to = n/2 — 5, where 5 = tan -1 (fr phot /r)). 
Finally, we assume that the flux and gas four-velocity profiles at 
the point of emission are given by their equatorial values F(r) 
and u'Xr), enabling us to calculate the spectra. Table 6 indi- 
cates that the effect on the estimates of a 4 may be significant at 
high spins. However, it is encouraging to note that the errors in 



the spin estimates decrease when we use a finite photospheric 
height. 

We should also note that Li et al. (2010) find a much larger 
effect on the spin estimates when they take the effect of a fi- 
nite photospheric height into account. We believe this is due to 
the fact that the \h/r\ profile in their analysis drops relatively 
sharply at small radii, like the profiles from the slim disc model 
and BHSPEC shown in Fig. 6. This disc geometry leads to self- 
shadowing of the disc. For our analysis in this subsection, on the 
other hand, we have chosen a constant \h/r\. Even if we were to 
use the \h/r \ profile from the GRMHD simulations, we would not 
expect to see significant self-shadowing, since the \h/r \ profile in 
the simulated model is nearly constant with radius (Fig. 6). 



4 A COMPARISON OF OBSERVATIONAL AND 
MODEL-DEPENDENT ERRORS 

The obvious question at this point is how big the errors in the 
spin estimates listed in Table 1 are compared with the observa- 
tional uncertainties in spin determination. We address this ques- 
tion in this section. The spins of eight stellar-mass black holes 
have been measured so far using the continuum-fitting method. 
The observational error estimates for the first four (see Shafee 
et al. 2006 for GRO J1655-40 and 4U 1543-^17, Davis et al. 2006 
for LMC X-3, and McClintock et al. 2006 for GRS 1915+105) are 
very approximate, and we disregard these results here. In more 
recent work, the principal sources of observational errors, as 
well as the uncertainties in the key model parameters (e.g., the 
viscosity parameter a), have been treated in detail. Moreover, in 
a recent paper on XTE J1550-564, Steiner et al. (2010) have ex- 
haustively explored many additional sources of error (see their 
Table 3 and Appendix A). The upshot of the work to date is 
that in every case the uncertainty in a, is completely dominated 
by the errors in three key dynamical parameters that are input 
when fitting the X-ray spectral data (McClintock et al. 2006). 
These parameters are the distance D, the black hole mass M, 
and the inclination of the inner disc i (which is assumed to be 
aligned with the orbital angular momentum vector of the binary; 
Li et al. 2009). In order to determine the error in a t due to the 
combined uncertainties in D, M and i, Monte Carlo simulations 
are performed assuming that these parameters are normally and 
independently distributed (e.g., Gou et al. 2009). 

Table 7 gives selected observational data for four black 
holes (all of these have been subjected to the rigorous error 
analysis described above): the inclination angle, which has an 
important effect on the model results (Tables 1-6); the spin pa- 
rameter; the absolute and fractional errors in r ISC0 (compare 
Table 2); and the luminosity. All errors are quoted at the 68% 
level of confidence. Note that the values of a„ range widely from 
~ to ~ 0.9. As a rough characterization, the uncertainties in 
the values of a, are Aa, ~ ±0.05 for the rapidly spinning pair 
of black holes and Aa t ~ ±0.2 for the slowly spinning pair. 
The corresponding fractional errors in r ISC0 range from approx- 
imately 10% to 20%. Comparing the fractional errors in r ISC0 in 
Table 7 with the closest counterpart results in Table 2 (i.e., clos- 
est matches for i and aj, we find that the error in the NT model 
is in all cases less than the observational error: A0620-00, 5.4% 
vs. 11.5%; XTE J1550-564, 8.2% vs. 17.9%; M33 X-7, 9.0% vs. 
10.7%; and LMC X-l, 1.9% vs. 20.5%. 

Furthermore, the estimates of the modeling error due to de- 
viations from the NT model obtained in this paper are very likely 
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Table 7. Data for four black holes: The entries are respectively the number of spectra analyzed; inclination angle; spin parameter; the approxi- 
mate/symmetrized absolute and fractional errors in r ISC0 ; and the Eddington-scaled luminosity. 



Black Hole 


No. 


i (deg) 


a* 


Ar(Ar/r) 




Reference 


A0620-00 
XTE J1550-564 
M33 X-7 
LMC X-l 


1 

45" 
15 
18 


51.0 ±0.9 
74.7 ±3.8 
74.6 ± 1.0 
36.4 ±2.0 


0.12 ±0.19 

34+ - 20 

u -0.28 

0.84 ±0.05 
92+ 05 


±0.65(11.5%) 
±0.86(17.9%) 
±0.29(10.7%) 
±0.43(20.5%) 


0.11 
0.05-0.30 
0.07-0.11 
0.15-0.17 


Gou et al. (2010) 
Steiner et al. (2010b) 
Liu et al. (2008, 2010) 
Gou et al. (2009) 



" Typical value: Number of spectra analyzed varies depending on details of data selection (see Steiner et al. 2010b). 



overestimates because the GRMHD simulation results necessar- 
ily correspond to relatively luminous discs: L/L Eii = 0.4 — 0.7, 
whereas the observed luminosities are typically only L/L Eii ~ 
0.15 (Table 7) and are strictly limited to L/L Edi < 0.30 (McClin- 
tock et al. 2006). Because the NT model improves as the thick- 
ness and luminosity of the disc decrease (Table 5), we conclude 
that use of the NT thin-disc model does not limit our accuracy. 
Rather, it is the uncertainties in the input parameters D, M and i 
that strongly dominate the error in a t . 



5 CONCLUSIONS AND DISCUSSION 

The main conclusion of this paper is that observational errors 
in current measurements of black hole spin by the continuum- 
fitting method dominate over the errors incurred by using the 
idealized NT model. We reached this conclusion by using 3D 
GRMHD simulations of thin discs to obtain realistic disc temper- 
ature profiles, then calculating the corresponding spectra, and 
finally fitting these spectra using the standard XSPEC model 
KERRBB. For disc thicknesses \h/r\ ~ 0.04 - 0.08, the errors 
in a, are up to about 0.2, depending on the inclination, for a 
non-spinning black hole, and up to about 0.1, 0.03 and 0.01 
for black holes with spins of 0.7, 0.9 and 0.98 respectively (Ta- 
ble 1). The errors in the spin estimates are particularly large at 
low spins and high inclinations, e.g., we find a spin estimate of 
0.15 for a non-spinning black hole viewed at an inclination angle 
of 75°. The results are quite close to those obtained by Reynolds 

6 Fabian (2008) for the iron line fitting method (see Fig. 5 of 
their paper) . Interestingly, we find that the fitted accretion rates 
are correct to within a few percent. 

A new and important contribution in this paper is that we 
establish via the slim-disc and BHSPEC models an approximate 
correspondence between the disc thickness as calculated from 
GRMHD simulations (§2.1) and the key disc observable L/L Eii 
(§3.1). Even though the simulated discs considered in this paper 
are geometrically quite thin, \h/r\ ~ 0.04 — 0.08, nevertheless 
it turns out that such discs correspond to fairly high luminosi- 
ties, I/I E dd ~ 0-4 — 0.7. For comparison, in observational work 
based on the continuum fitting method, the data-selection cri- 
terion L/L Eii < 0.3 (McClintock et al. 2006) is generally em- 
ployed (which, from Fig. 6, corresponds to \h/r\ < 0.03). The 
validity and usefulness of this criterion can be best judged by ex- 
amining the results for 411 observations of LMC X-3 in Steiner 
et al. (2010a, see their Figs. 2 and 3). For L/L EiA < 0.3, the 
inner disc radius r in is very nearly constant, rising only slightly 
at luminosities above L/L Eii ss 0.2. However, r in increases quite 
abruptly as the luminosity exceeds 30% of Eddington. Remark- 
ably, at these higher luminosities there is little scatter in the data 
and the increase in r in is smooth and systematic. 

It is difficult to say at this stage what the reason is for the 



above increase in the inner disc radius above the critical lumi- 
nosity of ~ 0.3L Edd . Both the GRMHD and the slim-disc mod- 
els predict that the inner disc radius should decrease (see the 
discussion of the radiation edge in Abramowicz et al. 2010) 8 . 
On the other hand, Li et al. (2010) were able to reproduce the 
observed increase by considering self-shadowing of the disc as 
a result of the off-midplane location of the disc photosphere. 
Interestingly, Abramowicz et al. (2010) find that the inner ra- 
dius of the disc is fairly close to the NT value for luminosities 
I < 0.3L Edd , and that the inner edge decreases quite abruptly at 
higher values of M. This, combined with the observed behaviour 
of r in in LMC X-3, may be a hint that something qualitatively dif- 
ferent happens at I/L Edd f=» 0.3; perhaps energy advection or 
disc self-shadowing becomes suddenly more relevant. 

One firm conclusion can be drawn from the results pre- 
sented in this paper. Since Table 5 indicates that the modeling 
error decreases as the disc thickness decreases, whatever the be- 
haviour of r in maybe above L/L Eii p» 0.3, at luminosities appro- 
priate to the continuum-fitting method (L < 0.3L Edd ), where the 
disc will be geometrically very thin, the errors in the spin esti- 
mates will be even smaller than those quoted in Table 1. There- 
fore, these errors are not a concern for the continuum-fitting 
method of measuring black hole spin. 

We must note one caveat about comparing the density scale 
height from the GRMHD simulations with that from the slim-disc 
and BHSPEC models: the latter models do not include magnetic 
pressure. The increase in the photosphere height due to mag- 
netic pressure could be as large as a factor of 2 (Hirose et al. 
2006), although other studies have found more modest changes 
(Blaes et al. 2006; Davis et al. 2009). 
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8 The caveat, as Abramowicz et al. (2010) point out, is that there are 
various ways of defining the inner edge of the accretion disc, and for 
some definitions, the inner disc radius can increase when the luminosity 
increases beyond ~ 0.3I Edd if the viscosity parameter is large enough 
(a > 0.2). However, the values of a that we see in our simulations are 
smaller, so this caveat does not present any problem. 
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APPENDIX A: MATCHING A GRMHD MODEL TO THE PAGE 
& THORNE SOLUTION 

Page & Thorne (1974, hereafter PT) define in their eqs (31a,b) 
two quantities: (i) /(r), which is proportional to the local disc 
flux F com (r) emitted from one side of the disc, as measured in 
the comoving frame of the fluid 9 , and (ii) w(r), which is pro- 
portional to the shear stress WT: 



f(r) = 4nrF com /M, 
w(r) = 2nrW*/M. 



(Al) 
(A2) 



These quantities enter the angular momentum and energy con- 
servation laws via (see PT, Eqs. 32a,b) 



(L f -w\ r =fL\ 
(E f -nw) r =fE f , 



(A3) 
(A4) 



where I'(r), E'[r) and fi(r) are the specific angular momen- 
tum, specific energy-at-infinity and the angular velocity. Eq. A3 
and Eq. A4 lead to the further relation (PT, Eq. 33) 



/ = -n r (£ t -ni t )- 



(A5) 



For a thin accretion disc, PT show that / has a general so- 
lution of the form 



/( r ) = -n r (£ t -^ 1 )" 



(E 1 -nL 1 )L + r dr + C 



(A6) 



9 The procedure for transforming the flux between the BL and comoving 
frames is described in §B2. 
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where C is an integration constant. They assume that the stress 
vanishes at the ISCO, and thereby determine the value of C. 
They then obtain the following particular solution: 



/ PT (r) = -n r (£ t -^ t )" 



(£* - m t )I t dr. 



(A7) 



Eq. 15n in PT gives an explicit analytical expression for / PT (r). 

We are interested in the following more general problem. 
We have a GRMHD numerical solution of a thin disc that has 
reached inflow equilibrium out to some radius r ie (defined in 
Penna et al. 2010). Beyond this radius, however, we cannot trust 
the numerical results, so we would like to match our simula- 
tion model to the PT solution. This will allow us to extrapolate 
the simulation beyond r ie and even beyond the radial range of 
the numerical grid. We wish to avoid the particular solution / PT 
given above since that assumes zero stress at the ISCO. Instead, 
we fit for the value of the integration constant C using the simu- 
lation. We also redefine the constant C slightly so that the fitting 
is made at r ie rather than r ISC0 . 

Let us write Eq. A6 as follows, 



/(r) = -n r (£ t -m t )" 



(£ f - m t )I t dr + C 



r > r ie 



(A8) 



where the new constant C is to be determined from the simula- 
tion at r = r iB . We can write 



(£ f - nL t )L t r dr 



(E f - nL f )L\dr 



(A9) 



(e 1 - ni 1 ) 2 



(e 1 - ni + ) 



/ptO) 



t\2 



(A10) 



/p T (>ie)- 



Substituting in Eq. A8, we obtain the result we seek: 

[{E f -ClL') 2 /n r ] r . 
f(r) = Mr) ~ ft^ ^ ' ^ /ptCO 



[(Et-nwn,], 



n, 



C, r>r i( 



where 



C = ■ 



(E + - fii + ) 2 

n~ r 



4nr ie F com (r ie )/M. 



(All) 



(A12) 



Except for F com (r ie )/M, the local disc flux of the simulation at 
r = r ie , all the other quantities are obtained from the idealized 
PT model and can be evaluated at any r outside the matching 
radius r ie . With this matching procedure, the local disc flux from 
the converged region of the simulation can be extrapolated to 
arbitrary radii. 

The choice of the matching radius r ie is an important issue. 
We would like to use as large a radius as possible, while still 
staying within the inflow equlibrium radius, but we need to en- 
sure that the choice of the matching radius does not affect the 
luminosity profiles d(L/M)/d(lnr) significantly. Fig. Al shows 
the luminosity profiles for a„ = 0, 0.7, 0.9 and 0.98 (bottom 
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Figure Al. Similar figure to Fig. 1, except that the solid line for each spin 
in Fig. 1 is replaced by a cluster of solid lines showing the luminosity 
profiles for different values of the matching radius, r ie /r ISC0 = 1.1, 1.2, 
1.3, 1.4 (solid lines, top to bottom). 



to top). For each spin, the cluster of solid lines shows the lumi- 
nosity profiles for r ie /r ISC0 = 1.1, 1.2, 1.3, 1.4, and the dashed 
line shows the NT profile. It is clear that for all the models ex- 
cept a, = 0.98, the choice of the matching radius has little ef- 
fect on the luminosity profiles. This turns out to be particularly 
true for r ie /r [sc0 < 1.3, as a closer look reveals. We therefore 
use r ie = 1.3r ISC0 in those three models. The a„ = 0.98 model, 
however, did not attain inflow equilibrium out to a sufficiently 
large radius. Its luminosity profiles therefore show more sensi- 
tivity to the matching radius. Hence we choose r ie = 1.2r ISC0 
for this model. This is a conservative choice, since it results in a 
larger deviation from the NT luminosity. The predicted errors in 
the spin estimates shown for this model in Table 1 are therefore 
likely to be overestimates. 



APPENDIX B: TRANSFORMATION BETWEEN THE 
BOYER-LINDQUIST AND COMOVING FRAMES 

Bl Transformation of vectors and one-forms 

To transform the four-momentum, we need to find the orthonor- 
mal basis of the comoving frame, which we do using the Gram- 
Schmidt orthonormalization procedure as described in Beckwith 
et al. (2008) and Shcherbakov & Huang (2011). Let us denote 
the comoving-frame basis vectors by e M and the fluid 4-velocity 
by u M . In the comoving frame, since the fluid is at rest, its four- 
velocity u is equal to e (t) . This is a coordinate-independent state- 
ment; it is true in any frame. We can thus denote the compo- 
nents of the e (r) in the Boyer-Lindquist (BL) frame by 

e H = (A 3 ,A 4 ,0, A 5 ), 



— (A 6 , A 7 ,A 8 ,A 9 ), 
e^ = (A 1; 0,0,^), 



(Bl) 



where the A; are determined by imposing orthonormality, i.e., 



(The index /! in is lowered using the BL 
metric g a g, while the index (v) is raised using the Minkowski 



metric 17 



diag(— 1, 1, 1, 1), since we want the comoving 



frame to be locally flat as well.) The e^ v) are the components of 
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the transformation matrix for four-vector components from the 
comoving frame into the BL frame. Some algebra gives 

e£) = (u\ u r , u e , u*), 

e[ r) = (u r u c , -(uy+u^u*), 0, u r u<f) I N r , 

e te) = (- u eU c , u e u r , l + u e u e , u e u*) / N e , 

e p ((j>) = (uj,, 0, 0, -u t ) I N^, 

and its inverse 



The total energy emitted in all directions is 
E = jdE = -eyjdp M 



(B9) 

coscS). (BIO) 



(-U t , ~U r , -Ug, -Us), 



e^> = (u r u c , -g rr (u t u t +u (^ ,li < ' , ), 0, u r it^) / N r , 

ef = -Asin 2 (u*. 0, 0, -u') / JV^, 
with 

N 2 = -g rr (u t u t +uXX 1+u e" 8 ) 

JV^-C^^+u^u^Asin^, 
A = r 2 + a 2 — 2Mr. 



The emission profiles we are interested in are isotropic (e^ = 
constant) or limb-darkened (e^ oc 2 + 3 cos 0) . For both of these, 
the second and fourth terms in the brackets vanish. The third 
term is multiplied by e[ e) oc u g = for the thin discs that we 
consider here, for which the fluid velocity is purely in the equa- 
torial plane. So we are left with 



= -u t \ 



d&en = -u t E a 



(BID 



(B2) 



Here, a is the (dimensionful) black hole spin: a = a t GM/c 2 . The 
transformation laws for vectors and one-forms are then given by 

^ = <%). ^) = ^v (B3) 
B2 Transformation of the flux 

To transform the radiation flux from the BL frame (which is what 
the GRMHD simulations give) into the comoving frame of the 
fluid, we use the following method. The comoving flux is defined 
as 

_ uiJ com "-^com 

u/1 com u L com v com 



Next, we need to transform the 3-volume d 3 V com . Con- 
sider first a four-volume in the BL frame bounded by dt, 
dr, d6 and d<p. The proper four-volume ij—gdtdrd6d(p = 
{^J—gf^dtdrdcpX^/gQgdO) is an invariant. (Here, g is the de- 
terminant of the covariant BL metric, and g cr ^ is the determi- 
nant of the t — r — cf> part of the metric.) Also, since the trans- 
formation into the comoving frame of the fluid involves a boost 
perpendicular to the 9 -direction, the proper length in the 9- 
direction (which is the term in the second set of brackets in the 
last expression) is also invariant. (This can be checked explicitly 
by examining the transformation of the vector = (0, 0, d 9, 0) 
into the comoving frame.) Therefore, the proper three-volume 
y/~Str<t>dtdrd<j> is also invariant. We therefore have 



d 3 V C om — V '-gtr<j>dtdrd4> = rdrd<pdt in the equatorial plane. 

(B12) 



(B4) 



The comoving flux then becomes 
dE 

Fcom ~ t-u t )rdrd(t>dt 
F 

-u, 



(B13) 
(B14) 



where d 3 V com is the 3-volume in the comoving frame. We first which is the required transformation for the flux, 
relate the energy in the fluid frame to that in the BL frame. Let 
ejj be the energy emitted per unit solid angle in the comoving 
frame. A ray of light emitted into a solid angle dQ. in a direction 
{9, <fi) (where 9 is measured with respect to the normal to the 
disc, i.e., the e m direction, and the c6 = direction is arbitrary) 
will then have an energy-momentum four-vector given by 



dpM = e^dCl(l,sm9 sin cJ, cos 9,sin9 cos c6), 



(B5) 



or, lowering the index, and remembering that the metric in the 
comoving frame is r) w(/3 j = diag(-l, 1, 1, 1), 



dp M = efjdn(-l,sin0sin<£,cos0,sin0cosc/>). 
Transforming this into the BL frame, we get 
d Pll = e^dp M . 



(B6) 



(B7) 



The energy-at-infinity of that ray as measured in the BL frame 
then is 



dE = -dp t = -e^dp (r) . 



(B8) 



APPENDIX C: RAY-TRACING GRID 

The grid in the image plane is generated using plane polar co- 
ordinates (b,/3). The radial grid points form a geometric series, 
b m = b q m , m = 0, . . . ,N b , where q is a number slightly larger 
than unity. This gives us high resolution close to the center of 
the image plane, which is necessary for resolving the inner re- 
gion of the accretion disc. The spacing in the angular direction 
is linear: /3„ = [3 + nSfi, n = 0, . . . , Ng. In addition, the grid is 
"squeezed" along the direction defined by the projection of the 
black hole spin axis onto the image plane (the "y"-axis), by an 
amount proportional to cos i, to account for the inclination of 
the observer. Therefore, the radial grid becomes a set of concen- 
tric ellipses; the Cartesian coordinates of the grid vertices are 
x mn = b m cosP n , y mn = b m sin/3„cosi. We choose i> = 1, which 
is inside the black hole shadow for all spins and inclinations, i.e., 
it is small enough that any rays with b < b fall into the black 
hole and can therefore be ignored. We also choose b Nb = 10000, 
which is large enough to correctly produce the low-energy end 
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of the spectrum in the energy range of interest (0.1 keV to 10 
keV). 

The flux in Eq. 2 then becomes 



1 

w _ 

cosi 

~w 

cosi 



I v dx dy 



I v bdbdp 



I v b 2 d(\o%b)dp. 



(CI) 



(C2) 



(C3) 



Since the values of log b and p are equally spaced, we can now 
use a higher order integration method like the Simpson's 1/3-rd 
rule to perform the integration (see, e.g., Abramowitz & Stegun 
1972). This gives considerably higher accuracy and faster con- 
vergence than a naive approach that simply sums up the product 
of I v in each grid cell with the area of that cell. 
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